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ABSTRACT 

The mechanisms and impact of correlated, or synchronous, firing among pairs 
and groups of neurons is under intense investigation throughout the nervous 
system. A ubiquitous circuit feature that can give rise to such correlations 
consists of overlapping, or common, inputs to pairs and populations of cells, 
leading to common spike train responses. Here, we use computational tools to 
study how the transfer of common input currents into common spike outputs 
is modulated by the physiology of the recipient cells. We focus on a key 
conductance — g^, for the A-type potassium current — which drives neurons 
between "Type II" excitability (low g^), and "Type I" excitability (high g^)- 
Regardless of g^, cells transform common input fluctuations into a ten- 
dency to spike nearly simultaneously. However, this process is more pro- 
nounced at low gA values, as previously predicted by reduced "phase" mod- 
els. Thus, for a given level of common input, Type II neurons produce spikes 
that are relatively more correlated over short time scales. Over long time 
scales, the trend reverses, with Type II neurons producing relatively less 
correlated spike trains. This is because these cells' increased tendency for 
simultaneous spiking is balanced by opposing tendencies at larger time lags. 
We demonstrate a novel implication for neural signal processing: downstream 
cells with long time constants are selectively driven by Type I cell populations 
upstream, and those with short time constants by Type II cell populations. 
Our results are established via high-throughput numerical simulations, and 
explained via the cells' filtering properties and nonlinear dynamics. 



INTRODUCTION 



Neurons throughout the nervous system — from the retina (phlens et al 



2008), thalamus (e.g., (Alonso et al., 1996)), and cortex (e.g., (Zohary et al. 
1994)) to motoneurons (Binder and Powers, 2001) — show temporal corre- 



lation between the discharge times of their spikes. This correlated spiking 



can impact sensory discrimination (Averbeck et al., 2006) and signal propa- 



gation (Salinas and Sejnowski, 2000). 



How do these correlations arise? We study a simple mechanism in which 
the inputs to a pair or population of neurons has a common component that 
is shared across multiple cells (Figure [l]). On an anatomical level, the large 
number of divergent connections that span layers and areas makes shared 



3 




Fig. 1: Shared input microcircuit, in which two neurons receive input cur- 
rents with a common component that represents correlated activity 
or shared afferents upstream. Each neuron is a single-compartment 
Connor-Stevens model (see Methods), with a maximal A-current con- 
ductance qa that we vary, eliciting a full range of Type I to Type II 
spiking dynamics. Shared input currents lead to correlated spikes, 
which are quantified as shown via spike counts ni, n 2 over sliding 
time windows of length T. The input currents received by each cell 
have mean [i and fluctuate with total variance a 2 ; the common noise is 
chosen with variance ca 2 , and independent noise terms with variance 
(1 - c)a 2 . 



afferents to pairs of nearby cells unavoidable (Shadlen and Newsome, 1998). 
Correlated spiking in areas upstream from the target cells can add to this 
anatomical factor. In fact, for some neural circuits, shared inputs are them- 
selves the dominant source of correlated spiking (Trong and Rieke, 2008). In 
general, correlating effects of shared input interact with effects of recurrent 
coupling (cf. (Ostojic et al., 2009)). 

What makes shared input circuitry especially interesting is the pivotal 
role of spike generating dynamics. For a given fraction of shared input, these 
dynamics control the fraction of shared output — that is, the fraction of 
spikes that will be shared across the two cells. This correlation transfer 
depends on two factors. The first is the mechanism of spike generation. The 
second is the operating point of the neurons (i.e., their rate and variability of 



firing (Binder and Powers, 2001 de la Rocha et al., 2007)). Excepting Hong 



and De Schutter (2008) , studies of correlation transfer have mostly focused 



on simplified neuron models, such as integrate-and-fire, phase, or threshold 
crossing systems, leaving open allied questions for models with more complex 
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subthreshold and after-spike dynamics. 

Here, we study correlation transfer for a family of conductance-based neu- 
ron models that varies from Type I excitability, in which firing can occur at 
arbitrarily low rates in response to a DC current (as for cortical pyramidal 
cells), to Type II excitability, in which firing occurs at a nonzero "onset" rate 



(as for fast-spiking interneurons or the Hodgkin-Huxley model) (Rinzel and 



Ermentrout, 1998; Hodgkin, 1948). We use the Connor- Stevens model (Con- 



nor and Stevens, 1971), which transitions between Type I and Type II as qa 

is varied 



sec 



- the maximal conductance of the A-type potassium current 
Fig. [2]). Beyond firing rates, Type I vs. Type II neurons differ in single-cell 
computation (Ermentrout et al., 2007) and synchronization under reciprocal 
coupling (Rinzel and Ermentrout, 1998). 

We test the hypothesis — based on predictions from simplified "normal 



form" phase models (Barreiro et al., 2010 Marella and Ermentrout, 2008 



Galan et al., 2007) — that the Type I to Type II transition will produce 
strong differences in levels and time scales of correlation transfer. Upon find- 
ing a positive result, and explaining it via the filtering properties of individual 
neurons, we ask how the distinct features of correlated spiking in Type I vs. 
Type II neurons manifest in signal transmission through feedforward neural 
circuits. Preliminary versions of some findings have appeared in abstract 



form (Barreiro, 2009 Shea-Brown et al., 2009). 



METHODS 



Circuit setup 

We primarily consider the feedforward circuit of Figure [TJ Here, each of two 
neurons receives two sources of fluctuating current: a common, or "shared" 
source I c = Gyfc^dt), and a "private" source I± = o\/\ — c£i(t) or I 2 = 



<T\/1 — c£ 2 (^)- Each of these inputs is chosen to be a scaled statistically 
independent, Gaussian white noise process (uncorrelated in time) — that is, 
(£i(£)£i(t + r)) = 5(t) for i = 1,2, c. This is for simplicity and agreement 



with prior studies of correlated spiking ( 


Lindner et al., 2005 


de la Rocha et 


al., 2007| 


Shea-Brown et al., 2008; 


Marella and Ermentrout, 2008 


Vilela and 


Lindner, 2009[ 


Barreiro et al., 2010 


). The common current J c has variance 



<j 2 c; each private current has zero mean and variance cx 2 (l — c). Note that 
these scalings are chosen so that the total variance of current injected into 
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each cell is always a 2 , while the parameter c gives the fraction of this variance 
that arises from common input sources. For example, when c = 0.5, 50% of 
each neuron's presynaptic inputs come from the shared and 50% from the 
independent input. Finally, the mean of the total current received by each 
cell is given by fi. This term represents the total bias toward negative or 
positive currents from all sources; in Fig. [TJ it is illustrated as part of the 
common input for simplicity. 
The combined currents, 

LppAt) = H + Ut) + Ii(t) 

= fj, + <ry/Z£ c {t) + aVl -c&{t) 

(i — 1,2) are injected into identical single-compartment, conductance-based 
membrane models (see Methods, "Neuron model"); spike times are identified 
from the resulting voltage trace. 



Neuron model 

We investigate correlation transfer in the Connor-Stevens model, which was 



designed to capture the low firing rates of a crab motor axon (Connor and 



Stevens, 1971 Connor et al., 1977). This model adds a transient potassium 
current, or A-current, to sodium and delayed-rectifier potassium currents of 
Hodgkin-Huxley type. The A- Type channel provides extended after-spike hy- 
perpolarization currents, which lead to arbitrarily low firing rates and hence 
Type I excitability (see Introduction and Methods, "Characterizing the dy- 
namics of spike generation," below). 
The voltage equation is 

C M ^ = -g L (V - E L ) - g Na m 3 h(V - E Na ) - g K n\V - E K ) (1) 



- g A A 3 B(V -E K )+I, 



app- 



The gating variables m, n, h, A, and B each evolve according to the standard 
voltage-gated kinetics; e.g., for m: 

dm m^iV) - m 

dt ~ r m (V) 1 ' 

where m^iV) is the steady-state value and r m (V) is the (voltage-dependent) 
time constant. All equations and parameters are exactly as specified as in 



I. 



H<Ibif M > I bif M<I b if M = I bif M > Ibif 




Fig. 2: Firing rate vs. injected current (/ — /) curves, for the deterministic 
(a = 0) Connor-Stevens model. Several values of Qa, yielding a range 
from Type II to Type I excitability, are shown — note the nonzero 
"onset" firing rates and Type II excitability for g& « mS/cm 2 , zero 
onset rate and Type I excitability for qa ~ 60 mS/cm 2 , and a grad- 
ual transition between. Insets show cartoons of dynamical transitions 
that lead to non-zero vs. zero onset rates: a subcritical Hopf bifurca- 
tion (left) and a saddle- node on invariant circle bifurcation (right). 



Connor et al. (1977) , with the exception that we vary the maximal A-current 



conductance, gA, over the range of values reported below. As qa is decreased 
from the value set by Connor et al. (1977) , the neuron transitions from Type 
I to Type II excitability; we describe this phenomenon next. 
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Measuring spike train correlation 

We represent the output spike trains as sequences of impulses yi{t) = 
ti), where t\ is the time of the kth spike of the ith neuron. The firing 
rate of the ith cell, (y^t)), is denoted To quantify correlation over a 
given time scale T, we compute the Pearson's correlation coefficient of spike 



counts over a time window of length T (as in, e.g., (Zohary et al., 1994 
Bair et al., 200l| )): 



Cov(ni,ri 2 ; 

Pt 



- V /Var(n 1 ) V /Var(n 2 ) 

where ni, n 2 are the numbers of spikes simultaneously output by neurons 1 
and 2 respectively in a time window of length T; i.e. rii(t) = f* +J yi(s) ds. 
If Pt is measured at values of T that are less than a typical inter-spike 
interval (ISI), we are essentially measuring the degree of synchrony between 
individual spikes. For larger T values, px assesses total correlation between 
numbers of spikes emitted by each cell. 



A short calculation (cf. (Bair et al., 2001 Cox and Lewis, 1966)) shows 
that Cov(ni,n 2 ) is 

Cov(ni,n 2 ) = T J C 12 (t)^lM dt (3) 

where the spike train cross-covariance Ci 2 (r) = (yi(t)y 2 (t + r)) — V\V 2 - Sim- 
ilarly, the variance Var(ni) can be given in terms of the spike train auto- 
covariance function. The autocovariance function of neuron 1, defined as 
Ai(t) = (yi(t)yi(t + r)) — uf, satisfies 

r T t — \t\ 
Var(ni) — T A x {t) 1 1 dt , 



— t T 



and similarly for neuron 2. 



Characterizing the dynamics of spike generation 

We give a brief discussion of spike-generation dynamics in the Connor-Stevens 
system. We organize this by focusing on how the neurons transition from 
quiescent (i.e., rest) behavior to periodic spiking in response to constant, 
noiseless (DC) current (i.e., when fluctuation terms o = 0). While these 
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results are well-established (Rush and Rinzel, 1995) we review critical aspects 
that will help us understand how correlated spiking arises in the circuit of 
Fig. [Tj 

For any value of g^, there is a critical value of the DC current, li = 
Iuf{gA)i such that the neuron has a stable stationary state — that is, it re- 
mains at a resting voltage — if /i < hif(gA)- If /i > hif(gA), then the neuron 
instead spikes periodically. We refer to these two regimes as subthreshold and 
superthreshold respectively. If the current fi is ramped slowly from a sub- 
threshold value, the membrane potential will shift slowly until /i = Ibif(gA), 
at which point periodic spikes begin; the rate continues to increase with /i 
over some range, as quantified by familiar firing rate-current (/ — I) curves 

(Fig.g. 

As discussed above, neurons are often classified as Type I vs. Type II 
based on whether their / — I curves are continuous vs. discontinuous (with 
a jump) at /i = I^f- Figure [2] shows that the Connor-Stevens model is Type 
II for gA ~ mS/cm 2 , Type I for « 60 mS/cm 2 , and displays a gradual 
transition in between. This qualitative shift in behavior is related to the un- 



derlying model (Equations [lpj) through bifurcation theory (Izhikevich, 2007 



Rinzel and Ermentrout, 1998[ ). The key fact is that all possible changes from 



resting to periodic spiking behavior fall into a few qualitatively equivalent 
categories of bifurcation. We next describe two of these categories that we 
will put to use below. 

In a saddle-node on invariant circle (SNIC) bifurcation, the stable resting 
state merges with an unstable resting state to create a single steady state 
precisely when \i = hif{gA)- A solution to Equations (TJ2) that begins near 
this steady state can return to it, via an excursion through state space that 
traces out a spike. This spiking trajectory persists for li > I bi f. However, 
due to the steady states at nearby values of /i, spiking trajectories take a 
long time to return to where they started — technically an infinite amount 
of time when fx = Iuf(gA)\ see the righthand schematic in Figure [2j As 
a result, the firing rate for \x near IufigA) is arbitrarily low, resulting in a 
continuous / — / curve and Type I excitability. 

In a subcritical Hopf bifurcation, by contrast, the stable resting state 
instead loses stability as an unstable periodic orbit shrinks into this point. 
For neural systems, there is typically also a stable periodic orbit in the state 
space. Once the resting state becomes unstable, trajectories quickly diverge 
from the resting state to the stable periodic orbit, which has a non-zero 
frequency fuf, see the lefthand schematic in Figure [2j This results in a 
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discontinuous / — / curve that jumps from zero (at rest) to fuf at lujiSA) 
— the characteristic of Type II excitability. 

We use software tools to automate the bifurcation analysis of the Connor- 



Stevens model: specifically XPP (Ermentrout, 2002), and MatCont (Dhooge 



et al., 2003). This allows us to track fixed points and limit cycles as system 



parameters g^, \i vary, and to check the mathematical conditions that define 
bifurcation types ( Guckenheimer and Holmes, 1983). 



Linking linear response theory, spike-triggered averages, 
and spike count correlations 

When the variance c of the shared input is small (see Fig. [I]), then we can treat 
the circuit with a shared input as a perturbation from two independently 
firing neurons. We describe this perturbation via linear response theory 



(Lindner et al., 2005; de la Rocha et al., 2007; Ostojic et al., 2009 Ostojic 



and Brunei, 2011), which is related to classical Linear-Poisson (LP) models 



of neural spiking (Perkel et al., 1967). That is, we make the assumption that 



the change in a neuron's instantaneous firing rate z/j(i) due to the shared input 
signal can be represented by linearly filtering the common (perturbing) input 



ui(t) = ( yi (t)\I c ) 

POO 

= f ,i + / K(s)I c (t - s) ds 
Jo 

= v Qji + (K*I c ){t) 



(4) 



where the filter K{t) = for t < (causality) and is it0 is the "background" 
average firing rate of the independently firing neuron. 

Eqn. ^ is extremely useful, because it isolates the common component of 
the response of neurons i — 1 and i — 2, which is an enhanced (or depressed) 
tendency to emit spikes, at a rate determined by the filtered, common input. 
As a result, the cross-covariance function, 



Ci 2 (t) = ((ui(t) -i/ ,i)Mt + r) -1/0,2)) 



becomes 
Ci 2 (r) 



OO /"OO 



K(s)K(s + r)(I c (t - s)I c (t -s + r-r))drds 
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This simplifies: because we assume J c is white (uncorrelated in time), (I c (t)I c (t+ 
t)) = Var(/ C )(5(r). Var(/ C ) is the instantaneous variance of the noisy process 
(here assumed to be constant), so that 



C 12 (r) = Var(/ C ) / K{s)K{s + r) ds 
Jo 

= Var(I c )(K*K)(T) (5) 



where K(t) = K(-t); cf. ( |Ostojic et al., 2009| ). 



Thus, the cross-covariance between our two spike trains has a simple 
expression in terms of the filter with which neurons process the common 
input into a firing rate. It remains only to identify this filter K(t). As in 



classical cases, this is precisely given by a spike triggered average (Gabbiani 



and Koch, 1998). Different from the classical setting, but as in Ostojic et 
al. (2009) , it is only the common component of the input current that is 
averaged in this procedure. 

The connection to the spike-triggered average may be seen by looking at 
the response to a white noise stimulus. On the one hand, 

(I e (t)ui(t + T)) = (I c (t)(K*I c )(t + r)} 

= (K*(I c (t)I c (t + r)))(r) 
= Vax(I c )K(r) 

if the noise is white. On the other, by ergodicity, averaging over repeated 
presentations and different stimuli I c will yield a response equivalent to aver- 
aging over repeated presentations of a single long (duration T max ) stimulus. 
As T max y oo, 



1 



Tmax 




(/ c (t)i/i(t + r)> = / I c (t)(y(t + r)\I c )dt 

h{t) (j2S(t + r-t k ) ) dt 

v 



,k=l 



T \ N 

- 1 max i 1 ' 



k=l 



(6) 



11 



which is exactly the average stimulus preceding a spike (multiplied by the 
firing rate). Therefore the linear response function K{r) is related to the 
spike-triggered average STA(r) as follows: 



K{r) - 



Var(/ C 
Var(/ C 



■STA(r) . (7) 



Below, we use this expression to derive K{j) from numerically computed 
spike triggered averages. 

Relating spike-triggered averages and spike-generating 
dynamics 

In order to relate the common input STA (defined in the previous subsection) 
to spike generating dynamics, it will be helpful for us to derive an explicit 
formula for the common input STA of a phase model, which captures the 
response of a tonically spiking neuron to a small-amplitude current I(t). Our 
formula, and the calculation that yields it, is very similar to a relationship 



previously derived (Ermentrout et al., 2007) for the STA of a phase oscillator 



without background noise. 

We consider a model which tracks only the phase of a neuron as it pro- 
gresses along its periodic spiking orbit: 

d9 

- = u + Z(9)I(t), 9e[0,2n) (8) 



The function Z(9), called a phase response curve or PRC ( 


Ermentrout and 


Kopell, 1984 Winfree, 2001 ; jErmentrout and Terman, 2010; 


Reyes and Fetz, 



1993), determines how a brief current injection applied at a specific phase of 
the cycle affects the timing of the next spike. By convention, the neuron is 
said to "spike" when 9 crosses 2n. 

To begin, we assume that the phase model is forced by scaled zero-mean, 
stationary stochastic processes, which we also label £ c (i) and For now, 

£ c (t) an d have unit variance and are differentiable with some finite 
correlation time r, although we will consider the limit r — > (i.e. the white 
noise limit). We are interested in the average value of £ c (i) that precedes 
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a spike; the term £j(t) will play the role of a background noise. Assuming 
that the background noise process is scaled by a small constant a, and that 
£ c is scaled by an order of magnitude e smaller still, we write the evolution 
equation of the phase model as 

d9 

— = l + Z(9)(a^(t)+aeUt)), 0e[O,T) 

where a, e <C 1 . Note that we have chosen our phase variable to have unit 
speed; i.e 9 G [0,T), where T is the period of the unperturbed (er = 0) 
oscillator. We proceed as in Ermentrout et al. (2007)] writing 9 as a series 
in the small parameters a and e 

9{t) = 9 (t) + a9 w (t) + e9 01 (t) + cr 2 9 20 (t) + ae9 11 (t) + e 2 9 02 (t) + ... 

and matching terms of same order in the evolution equation, we find 9 (t) = t. 
We additionally find that 6*01 = 0, 6*02 = 0, and 

9' 10 = Z(t)Ut) 

9' 20 = Z'(t)£i(t)9 10 (t) 

9' n = Z(t)Ut) 



so that 



ft 

9 w (t) = I Z(s)Ci(s)ds 



9 20 (t) = [ Z'(s)Us) f Z(r)Ur)drds 
Jo Jo 

9 u (t) = [ Z(s)Us)ds. 
Jo 

In order to compute the spike-triggered average, we need to find the time 
of the next spike, assuming the neuron has just spiked (0(0) = 0); in other 
words, the time r when 9{r) = T. As above, we expand 

t = T + otio + ct 2 t 20 + crer u + ... (9) 



Using our previous expressions for 9(r), and using the fact that r = T + 
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ar w + 0(a 2 , ae) to decompose the stochastic integrals, we find 

no = - / Z(s)£i(s)ds 
Jo 

T 20 = - / Z(s)Us)ds- Z'(s)Us) Z{r)Ur)drds 

Jt Jo Jo 

rii = - / Z(s)£ c (s) (is . 
./o 

Next, we use Taylor's theorem for smooth functions to expand £ c about T — t 
to compute 

STA(i) = (I c (r-t)) 

= ore(£ c (T + cttio + or 2 r 20 + crern - t) 

= ae(UT -t) + (<jt w + a 2 r 20 + <7ern)£(T - t) 

+\(vt 10 + a 2 r 20 + aemfCAT - t)) 

= ae(UT -t) + (<jt w + a 2 r 20 + <7ern)£(T - t) 

+ \° 2 rloC(T-t)) + 0{a\a 2 e). 

where we have kept terms up to second order both in our expression for r, 
and in our Taylor expansion of £ c . We can use the independence of £ c and & 
to eliminate a large number of terms, as 

Similarly, 

(um(t+ s )) = (ut))(&t+ 8 )) = o, (io) 

and so forth for expressions with higher derivatives. The only term that 
survives is 

STA(t) = ae(C(T - t) x -ae [ Z(s)£ c (s)ds) 

Jo 

= -M 2 f Z(s)(e c (T-t)Us))ds 
Jo 

= -{ae) 2 f Z{s)A' c {T -t- s)ds 
Jo 

= -(ae) 2 [ Z'(s)A c (T-t-s)ds, 
Jo 
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where A c is the autocovariance function of £ c and we used integration by parts 
in the final step. Taking the white noise limit (A C (T — t — s) — > 5(T — t — s)) 
and using the periodicity of the PRC {Z(T — t) — Z(—t)), we recover a very 



similar expression as Ermentrout et al. (2007) 



STA(t) 



ere 



'z'(-t). 



(11) 



Numerical simulations and estimates of spike count 
statistics 

To compute spike count correlations and other statistical quantities, we per- 
formed Monte Carlo simulations of the circuit in Fig. [T] The governing 
Connor-Stevens equations ([TJ [2]) were integrated using the stochastic Euler 
method with time step At = 0.01 ms, for a total time T max of 8 x 10 6 ms. Ran- 
dom input currents were chosen at each time step using a standard random 



number generator (Marsaglia and Zaman, 1994). To facilitate exploration 



of parameter space /i, a, c, and 5% we distributed computations on paral- 
lel machines through the NSF Teragrid program (htt p: / /www.teragrid.org] ). 
The simulation code was implemented in FORTRAN90, and distribution 
scripts in Python for running on clusters with and without PBS submission 
protocols. All code and scripts will be available at the modelDB site upon 



publication ( http: / / senselab.med.yale.edu / modeldb / ) . 



We register spikes in our simulations at times when the membrane voltage 
exceeds —30 mV and maintains a positive slope in voltage for the next three 
time steps (0.03 ms). To avoid counting each spike more than once, we omit 
a 2 ms refractory period after each spike. 

Spike count statistics were computed directly from the recorded spike 
times, based on a single long simulation, after discarding an initial transient 
(200 ms). When sampling spike counts over a time window T, we advance 
the window by \T, resulting in approximately 4T max /T (correlated) samples; 
consequently, our estimates of spike counts become noisier as T increases. 
To estimate standard errors on spike count statistics, we further divided 
the simulation into 10 equal time intervals (8 x 10 5 ms each) and computed 
statistics on each sub-simulation; the standard deviation, divided by VlO, 
gives us an estimated standard error of the mean. When appropriate, these 
are presented along with the mean estimates, as error bars. 

Below, we also report spike triggered averages (STAs) described above; 
these were computed using long simulations of length 8 x 10 7 ms for several 
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sets of parameter values fi, a, and g&. To compute these, the common 
input current I c was treated as the "signal" which was averaged and the 
private input as a "background" which was not. We used c = 0.10 in this 
computation. In our code, the history I c was continuously recorded for a 
duration into the past; when a spike was recorded, the STA was augmented 
by this current. 

Finally, we generate auto- and cross-correlograms (shown in Figure [4]) 
by collecting inter-spike intervals (ISIs) from our simulations in 1 ms-long 
bins. These are used, after the standard normalization, as auto- and cross- 
covariance functions. 



RESULTS 

Rich structure of spike count correlations over short and 
long time scales 

Our central findings contrast how different conductance-based neuron models 
produce correlated spiking when they receive overlapping fluctuating inputs, 
via the shared-input circuitry in Fig. [TJ Specifically, we show how this cor- 
relation depends on the Type I vs. Type II excitability class of a neuron 
described by the well-studied Connor-Stevens model. As discussed above, 
neurons are often classified as Type I vs. Type II based on whether their fir- 
ing rate-current curves are continuous vs. discontinuous at /i = J^/. Figure [2] 



demonstrates — as shown in (Rush and Rinzel, 1995) — that the Connor- 



Stevens model is Type II when the maximal A-current conductance qa ~ 
mS/cm 2 , Type I for ~ 60 mS/cm 2 , and displays a gradual transition in 
between. Thus, we fix the neurons in the shared-input circuit to a point 
along the spectrum from Type I to Type II excitability by choosing different 
values of g^. 

To compute levels of correlated spiking, we then fix the correlation in 
the input currents — that is, the fraction of the current variance that is 
shared vs. private to the two cells — to a preset value c. For each value of 
c and gA, we compute spike count correlations for wide range of operating 
points for the neurons, as determined by a ~ 200 x 50 grid of values for the 
mean current /i and variance a (both [i and a are sampled at a resolution of 
0.1 /iA/cm 2 ). Specifically, we vary fi over values centered at the threshold 
current Ibif(gA), from a minimum /i = Ibif(gA) — 10 (/iA/cm 2 ) to a maximum 
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M (|jA/cm 2 ) 

Fig. 3: Fano factor of spike counts over a long time window (T = 256 ms) for a 
~ 200 x 50 grid of values for the mean current \i and variance a. From 
top to bottom, Type II to Type I: (a) gA = 0, (b) gA = 30, and (c) 
gA = 60 mS/cm 2 . Markers indicate relative location of (fi, er)-pair; 
subthreshold by 1 /iA/cm 2 (diamond), superthreshold by 2 /zA/cm 2 
with low noise (circle) and high noise (square), and superthreshold 
with matched Fano factors (asterisk, see text). 

/i = Ibif(gA) + 10 (/xA/cm 2 ) for each value of gA- This enables us to cover, 
respectively, both subthreshold (i.e., fluctuation-driven, ji < I bi f(gA)) and 
superthreshold (i.e., mean-driven, /i > Iuf(gA)) firing regimes for each value 
of gA- We additionally vary a over < a < 5 (/zA/cm 2 ), so that we cover 
the range from nearly Poisson, irregular spiking to nearly periodic, oscillatory 
spiking. This is demonstrated by Figure |3j which shows the Fano factor of 
spike counts over a long time window (T = 256 ms) — a proxy for the squared 
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Fig. 4: 



Spike count correlations for three models at both short and long time 
scales. Each row displays data from a value of gA'- from top to bot- 
tom, Type II {qa = mS/cm 2 ), intermediate (gA = 30 mS/cm 2 ), and 
Type I (gA = 60 mS/cm 2 ). Left column: Spike count correlations 
PTi for short windows T = 4 ms. Center: Spike count correlations 
Pt, for long windows T = 128 ms. Markers indicate points used for 
cross-model comparison: subthreshold by 1 pA/cm 2 (diamonds), su- 
perthreshold by 2 pA/cm 2 and low noise (circles), superthreshold by 
2 pA/cm 2 and high noise (squares), and superthreshold with matched 
Fano factor (stars). Right: Cross-covariance and autocovariance (in- 
set) functions for the superthreshold high noise points (squares). Be- 
hind cross-covariance functions, the shape of the triangular kernel 
that relates this function to spike count covariance (as in Eqn. [3]) is 
illustrated for T = 4 ms (green) and T = 128 ms (yellow). For each 
value of gA, autocovariance functions are given in normalized units 
(so that A(0) = 1). 
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inter-spike interval coefficient of variation (Gabbiani and Koch, 1998) — over 
the entire p, a parameter space for three representative values of gA (0, 30, 
60 mS/cm 2 ). For each value of gA, the Fano factor spans a range from near 
zero (periodic) to one (i.e. Poisson-like) or higher. 

For each set of parameters gA, c, p and a, we compute the Pearson's 
correlation coefficient p? between the spike counts that the neuron pair in 
Fig. [T] produces in time windows of length T. Figure [4] summarizes the 
results, for inputs with 10% shared variance (c = 0.1). Here, we view over 
the entire p, a parameter space for three representative values of gA (0, 30, 
60 mS/cm 2 ) and two different time windows (T = 4 ms and T = 128 ms). 
Values of px depend in a strong but systematic way on all of the parameters 
we have introduced. As we move down a column, we see major qualitative 
differences in levels of correlation that emerge at different points through the 
Type II (gA = 0) to Type I (gA = 60) spectrum. Within each panel, the 
operating point set by input mean and variance (p, a) have a strong impact 
on pt- Finally, the levels and trends in pj- depend strongly on the time scale 
T. We now describe these trends in more detail; the sections that follow will 
give an explanation for how they arise. 

We begin with the upper panels in Fig. |4j which show correlation px 
for gA = and hence Type-II excitability. First, note that correlations are 
overall quite weak. The largest values of pt obtained are f» 0.04, indicating 
that ~ 40% or less of correlations in input currents are ever transferred into 
correlations in output spikes. Moreover, the level of correlations pt and their 
dependence on input parameters p and a appear roughly similar for both 
short and long time scales T. In both cases, for a fixed value of DC input 
p, a general trend is that pt gradually increases with fluctuation strength a. 
For a fixed value of a, in general pt first increases and then decreases with p; 
the dependence is slightly more complex at longer T. Significantly non-zero 
values of pr are present for p < I^f, as a becomes appreciably high; this 
reflects the bistable firing dynamics of the underlying deterministic system, 
which supports both a stable resting state and a stable spiking trajectory for 
H < hif (see Methods, "Characterizing the dynamics of spike generation"). 

For Type I excitability at g& = 60 (lower panels in Fig. [4]), the picture is 
dramatically different. First, there is a marked difference between correlation 
elicited on short vs. long time scales T, with much stronger correlations 
observed for larger T. Moreover, correlations produced by Type I neurons 
over longer time scales T are much higher than those observed for Type 
II neurons at any time scale: the largest values of pt obtained for Type 
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I are f» 0.8, indicating that ~ 80% of correlations in input currents can 
be transferred into spike correlations. Conversely, correlations for Type I 
neurons are strongly suppressed at short time scales, where ^ 10% of input 
correlations are transferred. Overall, trends in as p and a vary are similar 
to those found previously: correlations increase with a, and first increase - 
then decrease — with p. 

Correlation transfer in the intermediate model, qa = 30, displays trends 
between those of the Type I (gA = 60) and Type II (gA = 0) cases. As when 
gA = 60, spike count correlations pt are very low for short time windows T 
and attain intermediate to high values ~ 60% of input correlations trans- 
ferred for longer T. As for both = 60 and g& = 0, increases with noise 
magnitude a and displays a nonmonotonic trend with mean current p. 

We obtain additional insight into how spike count correlations depend 
on Type I vs. Type II spike generation and the time scale T by choosing 
matched values of input parameters p and a, and comparing spike count 
correlations produced for different values of the A-current conductance gA- 
We first concentrate on the p and a values indicated by squares and circles in 
Fig. |ij Both of these points indicate superthreshold inputs p = Iuf(gA) + 2 
(//A/cm 2 ) for all gA values. The square corresponds to higher noise a = 5 
pA/cm 2 , and the circle to lower noise a = 1 pA/cm 2 . In Figure [HJ we plot pt 
for a full range of T values from 1 to 200 ms, for nine values of g& between 
gA = and gA = 60 (thus filling in intermediate values of and T between 
those in Fig. For both superthreshold cases, we see that Type II neurons 
transfer more input correlation into output (spike) count correlation at small 
T, while Type I neurons transfer more at large T; this transition occurs, 
roughly, at a value T switch indicated by the dotted line. We note that for 
the low noise case a = 1 (Figure [51(b)), the trends appear less ordered as 
varies; as we will see in the next section, this is because the cross-covariance 
function is more oscillatory here, so that pr has not yet converged to its 
asymptotic large T value. 

Subthreshold points, denoted by diamonds in Fig. |4| were also compared: 
the mean input current is chosen to be p = I bi f(g A ) — 1 (/iA/chi 2 ), and the 
noise magnitude to be a = 5 pA/cm 2 . Here, the differing dynamical structure 
between Type II and Type I neurons is evident in the firing statistics (see 
Table [T]): while the bistable Type II neuron (gA = 0) sustains a substantial 
firing rate, the monostable Type I neuron (gA = 60) barely fires at this level of 
input current. The correlation coefficient pT, is also very low for g^ = 60 at all 
time windows (Figure pic)); this is consistent with the relationship between 
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Fig. 5: Correlation coefficient px vs. time window T. Colors indicate g& = 
(dark blue) through = 60 (red) mS/cm 2 . The data from the su- 
perthreshold cases (a, b, and d) show the switch from Type II cells 
transferring more correlations to Type I cells transferring more as 
T increases. The dotted line indicates the approximate time win- 
dow T switch where the switch occurs, (a) High noise, superthreshold 
(b) Low noise, superthreshold, (c) Subthreshold, and (d) high noise, 
where input current /i has been chosen to match spike train variability 
across different values of g& (see text). 



correlation and firing rate identified in earlier studies (de la Rocha et al., 2007 



Shea-Brown et al., 2008). Overall, note that the correlation coefficient px 



increases steadily with T for the Type I neurons (high g^) but stays roughly 
constant over a broad range of T for Type II neurons (low g^). Thus, while we 
do not observe a clear value of T switc h for all values of g& for the subthreshold 
point in Figure [5]^c), we see the same relative trends as for superthreshold 
points. Below, we will see how this effect follows from filtering properties of 
Type I vs. Type II cells. 
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Rpsimp 


Statistics 


Qa = 


04 = 30 


Qa = 60 


superthreshold, high a (□) 


i/ (Hz) 


113.2 
0.059 


69.9 
0.0795 


31.6 
0.195 


superthreshold, low er (0) 


v (Hz) 

FFt=256 


108.4 
0.0107 


65.0 
0.013 


34.8 
0.023 


subthreshold (0) 


v (Hz) 

FFt=256 


81.7 
0.145 


30.0 
0.282 


0.171 
0.99 


superthreshold, 
fixed variability (^K) 


v (Hz) 

FFt=256 


113.2 
0.059 


82.5 
0.059 


68.5 
0.059 



Tab. 1: Output firing statistics at each of the comparison points identified in 
Figure [4] (see text). For each set of matched points, we note the firing 
rate v and the Fano factor of spike counts over long time windows 
(specifically, T = 256 ms). 



Because spike generation mechanisms vary widely as changes, the neu- 
ron models with matched input statistics at different values of gA in Figs. [5|a, 
b, c) do not all have the same firing variability. In the final panel of Fig. |5j 
we address this by showing that the same trends in px persist if we select val- 
ues of p to maintain constant firing variability for each value of g A (see Table 
[T] variability measured via large-time (T = 256 ms) Fano factor). Here, we 
fix a = 5 pA/cm 2 ; the required current value p for each g^ is indicated with 
a red star in Fig. |4} 

In sum, for matched values of the mean and variance of input currents, 
a pair of superthreshold Type II (vs. Type I) neurons will produce greater 
spike count correlations pt at short time scales T. For a wide range of 
choices for the mean and variance, there will be a value of T switc h where this 
relationship reverses, so that Type I (vs. Type II) neurons produce greater 
Pt for T > T switc h- For matched subthreshold currents, similar trends are 
present; overall, the presence of a time T switc h depends on how the input 
statistics are chosen. 

Finally, we note that the general trends observed here carry over, largely 
unchanged, to different values of c. Figure [6] shows that, for the range 0.1 < 
c < 0.5, trends in how pt changes with excitability type via (gu) remain 
consistent. In particular the relationship between input correlation c and 
spike count correlation pj- is roughly linear over this broad range of shared 
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Fig. 6: Output correlation coefficient pt vs. input correlation coefficient c, 
showing an approximately linear relationship. Left: short time win- 
dow (T = 4 ms). Right: long time window (T = 150 ms). Colors 
indicate gA = (dark blue), g A = 30 (light blue), and g A = 60 
(red). Markers indicate relative location of (/x, <r)-pair; subthreshold 
(diamond), superthreshold with low noise (circle), and superthreshold 
with high noise (square). 



inputs. 

Trends in cross-correlation functions for Type I vs. Type II 
neurons 

The trends in spike count correlations that we have just described can be 
explained from the cross-covariance functions for neuron pairs, and how they 
differ as the characteristics of input currents and the level of the A-current 
conductance vary. We now demonstrate this via the cross-covariance 
functions shown in the right-most column of Fig. |4j these are for the su- 
perthreshold high noise cases (a = 5) discussed in the previous section. 

To make the connection, recall that the spike count covariance, Cov(ni, 7x2), 
measured over a window of duration T is given by the integral of the cross- 
covariance function Ci_2(t) against a triangular kernel of width T (see Meth- 
ods, and Equation [3]). Thus, for short windows T, only the central peak 
of C\2{t) contributes to spike count covariance. In the limit of long win- 
dows T — > 00, the spike count covariance is simply the integral of the cross- 
covariance function, multiplied by T, over the whole r axis. Spike-count cor- 
relation p T is then given by the ratio of spike count covariance to the spike 
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count variance. As we will show below, and ov ^, 1 ' n2j often show similar 
trends with T. Both quantities are of interest: while px gives a normalized 
metric of correlation, Cov (^ 1 ' n2 ) j s the relevant quantity to analyze impact 
on downstream excitable cells (see below, "Readout of correlated spiking by 
downstream cells"). 

Armed with these relationships between Ci 2 (r), Cov(n 1 , n 2 ), and px, we 
revisit the trends observed in the previous section and explore their origin. 
Starting in the upper-right of Fig. [5J note that C\2(r) has a much larger 
central peak — and hence short-T spike count covariance — for Type II 
excitability = 0) than for Type I (bottom-right, qa = 60). This is also 
clear in Fig. [7^ a), where we plot spike count covariance vs. T. Over long 
windows T, the trend reverses. For qa = 0, the cross-covariance function 
shows oscillations with significant negative and positive lobes. These lobes 
tend to cancel as C\2{t) is integrated over long windows T. This cancellation 
results in little overall change in values of spike-count covariance computed 
at increasingly long values of T. For Type I excitability, however, Ci2(t) is 
mostly positive, so that spike count covariance increases with T. 

As discussed above, spike-count correlation px is given by the ratio of 
spike-count covariance and variance. Comparing Fig. [7[a) and Fig. [5]^a) it 
is clear that spike-count correlation and spike-count covariance display the 
same trends as qa is varied. For example, for very short times T spike-count 
correlation px is given by the ratio of the peak in C\i(j) to that in A\{r)\ 
this ratio is also larger for Type II vs. Type I excitability. 

For other operating points (p, a), while trends in spike-count correlation 
and spike-count covariance do not exactly agree, the some relative trends 
persist. Specifically, the general pattern that Type II cells produce greater 
covariance over short time windows, and that this trend disappears or re- 
verses for larger time windows, holds for each of the superthreshold operating 
points explored here. Overall, the major trends in spike count correlations 
stem from the presence vs. absence of large negative lobes in cross-covariance 
functions for Type II vs. Type I neurons. We next describe how this differ- 
ence arises via the distinct filtering properties of the two neuron types. 
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Fig. 7: Spike count covariances, and their relationship to spike triggered av- 
erages (STAs). Each row compares data collected at a comparison 
point for the input current statistics (/i, a); see text. From top, su- 
perthreshold current with high noise (square), superthreshold with 
low noise (circle), subthreshold (diamond), and superthreshold with 
high noise and matched variability (red star). (Left column) Actual 
(dotted lines) and predicted (heavy solid lines) spike count covari- 
ances (Cov(rii, n 2 )/T), for representative points and all qa values. 
Colors identify qa values, which range from gA = (dark blue), 
through g A = 30 (light blue), to g A = 60 (red) mS/cm 2 . (Right 
column) Select spike-triggered averages (right column) and one-sided 
cross-covariance functions (left column, derived from the STA using 
Eqn. IB"]) used to compute predicted spike count covariances. 
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Spike triggered average currents reliably predict spike 
count covariance 

The previous section showed how trends in spike count covariances for Type 
I vs. Type II neurons follow from the presence of both strongly negative 
and positive lobes in cross-covariance functions for Type II neurons. Here, 
we explain the origin of this phenomenon. Equations ((5J7J) (see Methods) 
provide the key link, in which the cross-covariance function is given in terms 
of a cell's spike-triggered average (STA), which is an estimate of the filter 
through which cells turn incoming currents into time-dependent spiking rates. 
Here, we define the STA in response to common input as an average of the 
currents that precede spikes over a single long realization: 



STA = (^I>^- r )J' ( 12 ) 

where the tk are the N spike times from the realization. 

We first show that the prediction of spike count covariances from STAs is 
accurate. The left column Fig. [7] shows close agreement between spike count 
covariances computed from "full" numerical simulation (heavy solid lines) 
vs. predictions from STAs via Eqn. J5|) (dotted lines). Next, we examine the 
shape of the STAs themselves (Fig. Yu right column). 

For each operating point we consider, the Type II STA (g A = 0) has 
a pronounced negative lobe. Functionally, this corresponds to a "differen- 
tiating" mode through which inputs are processed: negative currents suffi- 
ciently far in the past tend to drive more vigorous spiking. Biophysically, 
this corresponds to the kinetics of ionic currents, such as inward currents 
that can be de-inactivated through hyperpolarization. By contrast, Type 
I STAs show a less prominent negative lobe, or none at all. The resulting 
filtering of inputs is characterized as "integrating:" a purely positive filter is 



applied to past inputs to determine firing rate (cf. (Dayan and Abbott, 2001 



Agiiera y Areas et al., 2003; Mato and Samengo, 2008)). 



The consequences for spike cross-covariance functions are straightforward. 
Trends are most pronounced for the superthreshold, high a case (Figure |7[a)). 
Here, the pronounced negative lobe in the Type II STA (g^ = 0) leads to 
a similar negative lobe in cross-covariance, and hence a sharp decrease - 
following an initial increase — of spike count covariance as a function of time 
window T. The Type I STA (g^ = 60) is positive, leading to a spike count 
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covariance which steadily increases until it overtakes the Type II value at 
T w 20 ms. These trends are also reflected in the spike count correlation px, 
as described previously. Moreover, analogous plots for superthreshold, high 
a points, where spike-count variability (Fano factor) is maintained across g& 
values (Fig. |7](d)), show the same trends. 

In the low a case (Fig. u£b)), the STAs have similar characteristics; how- 
ever there is signifiant ringing in the STA at the characteristic frequency of 
the oscillator. By the time the predicted (and actual) covariances |;Cov(ni, 77-2) 
reach a limiting value, they are close to zero, and possibly too variable to 
order definitively. It appears that covariance is still larger for Type II than 
for Type I at T = 200 ms. In the subthreshold case (Fig. [7]^c)), the STA for 
the Type I neurons is very small, consistent with the very low firing rate here. 
The Type II neuron shows a more robust response, similar in magnitude but 
less oscillatory than for the superthreshold regime. 

Trends in spike-generating dynamics mirror trends in 
spike-triggered averages and transferred correlations 

The transition from Type II to Type I spike generation in the Connor-Stevens 
model — as manifest in the progression from discontinuous to continuous 
spike- frequency vs. current curves in Fig. [2] — can also be characterized via 
the type of bifurcation that governs the transition from quiescence to periodic 
spiking as increasingly strong currents are injected (see Methods). 

For < gA < 46 mS/cm 2 , the transition occurs via a subcritical Hopf 
bifurcation, as voltage trajectories jump from a stable rest state to a pre- 
existing stable periodic orbit (limit cycle). This transition is schematized 
in the upper-left cartoon in Fig. [2j As this figure shows, for smaller values 
of gA in this range, the frequency of this cycle is high (> 60 Hz). The 
voltage-conductance dynamics near both stable structures - the stable rest 
state and the limit cycle - is oscillatory. This creates a resonator property 



(see (Izhikevich, 2007) and references therein): if they are properly timed, 
both negative and positive inputs cooperatively produce spikes or cause them 
to occur earlier than they would in the absence of inputs. This is reflected 
in the negative and positive lobes in the STA oc K{r) for the gA = cases 
in Figure [7j recall that the STA is the filter applied to incoming currents to 
determine firing rates. 

By contrast, for large gA > 58 mS/cm 2 , a saddle-node on invariant circle 
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(SNIC) bifurcation occurs (see Methods). As sketched in the upper- right 
cartoon in Fig. [2j in this case there is a pair of fixed points that form a 
"barrier" to spike generation for subthreshold values of fi, and the shadow, 
or "ghost" of these fixed points still affects dynamics for superthreshold \x 
- producing slow dynamics in their vicinity. The inputs that will elicit or 
accelerate spikes are those that will push trajectories past the fixed points, 
or their ghost, in a distinguished direction. These inputs therefore tend to 
have a single (positive) sign. This is referred to as an integrator property, 
and gives rise to the mostly positive STAs seen in Fig. [7] for the g& = 60 
cases. (This argument breaks down for very low-variance (low a) inputs, as 
we will see in the next section). 

Between these two extremes in the minimum frequency in response 
to a ramp current decreases steadily, creating a gradual shift between Type I 
and Type II behavior. This gradual shift is mirrored in the neural dynamics, 
in which the slow regions in the state space become increasingly dominant. 
This transition is clear in the spike triggered averages — and therefore spike 
count covariances — shown in Fig. [7j For example, for g^ = 30, we find 
both distinctly "Type II" -like and "Type I" -like aspects in the high noise 
(Fig. [7[ a)) and subthreshold (Fig. [7^c)) covariance trends respectively. In 
the former, spike count covariance increases — then decreases — with T; 
reflecting an oscillating STA; the end result is that the (normalized) covari- 
ance at T = 200 ms is lower than the covariance at T = 1 ms. In the latter, 
the normalized covariance steadily increases with T, reflecting a non-negative 
STA. 



Phase response curves (PRCs) predict common-input STAs 

We next show how, for superthreshold operating points, the key properties of 
Type I vs. Type II spike generation that determine the filtering of common 
inputs can be understood via a commonly used and analytically tractable 
reduced model for tonically spiking neurons. The response of such neurons 
to an additional small- amplitude current I(t) can be described by a phase 
model, a one-dimensional description which keeps track only of the progress 
of neuron along its periodic spiking orbit (or limit cycle). Identifying progress 
along the cycle with a phase 6 G [0, 2n), this model is completely determined 



by a single function of phase Z(8), called a phase response curve or PRC (Er 



mentrout and Kopell, 1984 |Winfree, 2001 Ermentrout and Terman, 2010| 
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Fig. 8: Comparison of PRCs to spike triggered averages (STAs) computed 
for both low and high noise, for both Type I (left column) and Type 
II (right column) neurons. For simplicity of visualization, each curve 
has been normalized by its maximum; that is Z(t) = Z(t)/ max(Z(i)), 
STA(t) = STA(t)/max(STA(t)), and ^Z'{t) = -Z'(t)/ max(-Z'(t)). 
In addition, the time axis has been scaled by the mean period in each 
case. Top row: PRC, showing monophasic and biphasic shape for 
Type I and Type II neurons respectively Middle row: High noise 
STA; the Type I neuron has lost the negative lobe in its STA, while 
the Type II neuron retains a negative component. Bottom row: Com- 
parison of (blue dashed) dPRCs with STA for the low noise (cyan) 
case. Both STAs have negative components. 
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Fig. 9: Correlation coefficient p at time window T = 200 ms, as Qa is varied. 
The data (gray solid) are from high-noise, superthreshold points and 
are the same as reported in Figure [7j The prediction (black solid with 
diamonds) uses Equation 14 These data show the increase in long 



time scale correlation as the model transitions from Type II to Type 
I. 



Reyes and Fetz, 1993): 



- = cu + Z(6)I(t). (13) 

We can interpret the meaning of this function by considering its effects on the 
timing of the next spike delivered at a particular phase of the limit cycle 0. If 
Z(4>) > 0, then a positive input delivered at that particular phase will push 
the neuron further along, advancing the time of the next spike; if Z(<f>) < 0, 
the same input would delay the time of the next spike. 



Neurons that display Type I spiking have a purely positive (or Type I) 
PRC, while Type II neurons show a PRC that has both positive and negative 



lobes ( 


Ermentrout and Kopell, 1984 


Ermentrout, 1996 


Hansel et al., 1995 


Brown et al., 2004). A purely posil 


ive PRC is characteristic of dynamics 



near a saddle node bifurcation, in which the system lingers near the ghost 
of its fixed points (as described in the previous section); input in a specific 
direction is needed to force the system away and elicit a spike. A biphasic 
PRC reflects oscillatory structure in the phase space, in which correctly timed 
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negative and positive inputs can cooperate to elicit a spike (as with a Hopf 
bifurcation). 

Strong relationships between the PRC and the STA have been found 
for neurons close to the threshold for periodic spiking (i.e., /i <^ I bi f, see 
Methods). Spike-triggered covariance analysis of both a Type I phase model 
and the Wang-Busaki model show that the dominant linear "feature" (corre- 



sponding to the STA) qualitatively resembles the PRC (Mato and Samengo 



2008) in the presence of sufficient current noise. In the (Type II) Hodgkin- 
Huxley model, the two dominant "spike-associated" features identified through 
covariance analysis closely resemble the STA and its derivative; the STA, in 



turn, closely resembles the PRC (Agiiera y Areas et al., 2003). 



In contrast, phase models in the oscillatory regime (far from the excitabil- 
ity threshold) are known to have an STA proportional to the derivative of 
the PRC (Ermentrout et al., 2007). We generalize this to the case of Fig. [Tj 



where the relevant signal £ c is delivered on top of a noisy background (see 



Eqn. gTJ, Methods): STA(t) oc -Z'(-t). 

In Figure |8j we test the accuracy of these relationships for the superthresh- 
old points considered above. We show results for Type I (qa = 60, left panels) 
and Type II neurons (qa = 0, left panels), and compare the STA computed 
at two different noise levels to the shape of the PRC (Z(8)) and its derivative, 
labeled dPRC (Z'(8)). The time argument of the STA has been scaled so 
that one period (T) maps onto the unit interval; likewise, the PRC is mapped 
onto the unit interval. At the lower level of noise, we have good correspon- 
dence between the STA and the dPRC in both cases. Notably, both Type I 
and Type II neurons have biphasic STAs. At high noise levels, while there is 
not a strong quantitative relationship between the STA and the PRC itself 



(unlike in the excitable regime explored by (Agiiera y Areas et al., 2003| )), 



the PRC carries important clues about the qualitative behavior of the STA. 
The Type II neuron retains the biphasic shape reflective of its PRC, while 
the Type I neuron has shifted to a purely positive STA. In sum, by predicting 
the STA shape, the PRC gives important clues to the linear response (and 
hence common input transfer) that we obse rve in Figure [7[ 



Finally, we test an alternate result (cf. (Barreiro et al., 2010 Abouzeid 



and Ermentrout, 2011)) that, in limited cases, relates PRCs to spike count 



correlations directly. For for long T and reasonably small a and c, 
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In Figure |9j we show that this gives a close approximation to simulation 
results for T = 200 ms in the superthreshold, high-noise case (see Table 1). 
Moreover, we can gain insight into the limitations this asymptotic approxi- 
mation by comparing with the superthreshold, low-noise case. The results 
of Barreiro et al. (2010) Abouzeid and Ermentrout (2011)] are derived by 
taking the asymptotic limit T — > oo before considering a finite but small; 
in practice, the smaller the noise variance a, the longer T must be in order 
to see this effect. For our low- noise points, the asymptotic behavior has not 
been recovered even at T = 1000 ms (as may be seen in Figure (5^b)). By 
using a very large (but probably biologically irrelevant) time window (data 
not shown), we eventually recover results consistent with the asymptotic 
prediction (Equation 14). 



Readout of correlated spiking by downstream cells 

How could the difference in spike count correlation between Type I and Type 
II cells impact neural circuits? We explore this impact in a simple network, 
in which correlated Type I or Type II cells collectively converge to drive a 
neuron downstream (see Fig. [T0[a)). 

In more detail, the drive comes from a population of N = 200 identical 
Type I (g A = 60 mS/cm 2 ) or Type II (g A = mS/cm 2 ) upstream neurons; we 
refer to these as population I and population II respectively. The upstream 
populations receive correlated inputs with c = 0.5 and values of // and a 
that yield matched levels of variability, as for the parameter set identified 
with stars in Fig. [4] (for population I, /i — 18 /xA/cm 2 and a = 5 /iA/cm 2 ; 
population II, [i = —6 /iA/cm 2 and o = 5 /iA/cm 2 ). This yields firing rates 
of uj = 63.5 Hz for neurons in population I, and ujj = 113 Hz in popula- 
tion II. Each upstream neuron has a single, instantaneous (delta function) 
synapse onto the downstream neuron of strength gj or gu; the relative size 
of the EPSPs are chosen so that the mean driving current is equal for each 
population [yigi = vngu, so that gi = 0.825 mV, gu = 0.5 mV). 

The total input received by the downstream neuron, I ds , is thus the 
weighted sum of A" upstream spike trains Uj{t): 

Id S (t) = gi^yj(t) or I ds (t) = gii^2yj(t) . (15) 

j j 

When the population size A^ is large, the summed signal has the same tempo- 
ral characteristics as the the cross-covariance between neuron pairs. Specifi- 
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tLIF< mS ) 1 L IF< m S) 



Fig. 10: (a) Schematic of "upstream" Type I or Type II neuron population re- 
ceiving common and independent inputs, and converging to a Leaky 
Integrate and Fire (LIF) cell downstream, (b) PSTHs from Type I 
and Type II upstream populations, (c) Predicted power of the volt- 
age fluctuations in the LIF cell, using STA (see text), (d) Actual 
firing rates of the LIF cell, showing similar trends with LIF time 
scale r LIF . 
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cally, the autocovariance of the summed input is 

A ds {r) = E[(I ds {t)-(I ds ))(I ds (t + r)-(I ds ))} 

= N(N - 1) E [( Vi (t) - Vi){ yj {t + r) - Vj )\ + NE[( yi (t) - H 

« iV 2 E [(!/<(<) + T)-!/.,-)] 

= iV 2 C 12 (r). 

This relationship is evident in the peri-stimulus time histograms (PSTHs) in 
Fig. [l~0|^b). For population II, fast fluctuations above and below the mean 
population output reflect the negative lobe in C\2{t) adjacent to its large 
peak. Meanwhile, fluctuations in the output of population I are less extreme 
and more gradual in time. 

The downstream cell integrates I ds (t) via Leaky Integrate-and-Fire (LIF) 



voltage dynamics (Dayan and Abbott, 2001): 



TLIF d ^ = -(V-V r )+I ds (t) 

where tlif is the membrane time constant and V r = —60 mV is the rest 
voltage. Spikes are produced when the voltage V crosses V t hresh = ~ 45 mV, 
at which point V is reset to V r . 

For the parameters we have chosen, the downstream neuron is driven sub- 
threshold, so that (I ds (t)) is not sufficient to excite a spike — any spikes must 
be driven by fluctuations in I ds (t). Thus, the variance of fluctuations in V(t) 
should give a rough estimate of how often membrane voltage will exceed the 
threshold, and consequently the downstream firing rate. This variance is easy 
to compute for a passive membrane (i.e., neglecting spike-reset dynamics). 
First, note that 



V(t) = V r + I I ds (s)L(t-s)ds 

J — oo 



where L is a one-sided exponential filter 



L(t) = -exp(-t/r L/F ), t>0 
r 

= 0, t < 0. 



We compute the variance as follows, using the causality of L to take each 
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upper limit of integration to infinity: 



E[V(tf 



E 



where L{t) = L(—t); the last step involved the substitution z = s — r and 
switching the order of integration. This final interior integral can be evalu- 
ated in the Fourier domain: using the properties that J r [L](a;) = F[L]{— u) 
and the fact that for real functions F[f](—oj) = J-[f)(uj), we find 
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(L*L) (to) 



\T[L\{u)) 
1 



1 + uHh F • 

Therefore (for example by consulting a transform table) 

1 



(L*L)(t) 



2t l 



exp(-\t\/T LI F) 



IF 



Substituting into Equation (18), we see that the variance of the downstream 



cell's voltage is given by a formula similar to that for the spike count covari- 
ances (Equation [3]): both involve integrating the cross-covariance function 
against a (roughly) triangular-shaped kernel, with time scale tlif in the 
former case and T in the latter. 

Figure |j~0[c) shows the by-now familiar trends that this predicts. For 
short membrane time scales t lif , Type II populations drive greater voltage 
variance; this is precisely analogous to the finding that spike-count correla- 
tions are greater for Type II cells over short time scales T. For long tlif, 
Type I populations drive greater voltage variance, just as Type I spike trains 
are more correlated over long time scales T. In panel (d), we compare this 
trend with actual firing rates elicited in the downstream cell (from numerical 
simulation). The general trends match, validating our simple prediction. 



(I ds (s) - (I ds })L(t -s)ds (I ds (r) - (I ds ))L(t - r) dr 

oo J — OO 

/OO POO 
/ dsdrL(t- s)L(t-r)~E[I(s)I(r)} 
-oo J — oo 

/OO POO 
/ dsdrL(t- s)L(t-r)A ds (s-r) 
-oo J —oo 

dzA ds (z)(L*L)(-z) (18) 
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In sum, downstream neurons with short membrane time scales (tlif ^ 5 
ms) are preferentially driven by Type II cells upstream; for longer time scales, 
the preference shifts to Type I cells. Some implications of this finding are 
noted in the Discussion. 



DISCUSSION 

Diverging connections, leading to overlapped input shared across multiple 
neurons, are a ubiquitous feature of neural anatomy. We study the interplay 
between this connectivity pattern and basic properties of spike generation in 
creating collective spiking across multiple neurons. We range spike generation 



over the fundamental categories of Type I to Type II excitability (Rinzel and 



Ermentrout, 1998; Hodgkin, 1948). The transition in excitability is produced 
by varying the A-current conductance within the well-studied Connor- 
Stevens neuron model. 

Our principal finding is that excitability type plays a major role in how 
shared — i.e., correlated — input currents are transformed into correlated 
output spikes. Moreover, these differences depend strongly on the time scale 
T over which correlations are assessed. At short time scales T, Type-II 
neurons tend to produce relatively stronger spike correlations for comparable 



input currents (Marella and Ermentrout, 2008 Galan et al., 2007). At longer 



time scales, the opposite is generally true: for a broad range of input currents, 
Type-I neurons transfer most of the shared variance in their inputs (~ 80%) 
into shared variance in output spikes, while Type-II neurons transfer less 
than half (~ 40%). 

We show that these results have direct implications for how downstream 
neurons with different membrane time constants will respond to Type I vs. 
Type II populations. Specifically, downstream neurons preferentially respond 
to populations that are strongly correlated on time scales similar to their 
membrane time constant. Interestingly, for the case we study, we find that 
the breakpoint between selectivity to Type I vs. Type II populations was for 
downstream membrane time constants of ~ 5 ms, easily within the ranges 
found experimentally. 

This raises interesting possibilities for neuromodulation. The membrane 
time constant of the downstream cell could be changed by shunting affects 
of additional background inputs — leading to a switch in its sensitivity to 
different upstream populations. Alternatively, modulators applied to the 
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upstream populations themselves could change their excitability from Type 
I to Type II (Stiefel et al., 2008 Stiefel et al., 2009), adjusting their impact 



on a downstream cell with a fixed membrane time constant. 

Overall, we demonstrate and apply a general principle: the presence and 
balance among different membrane currents controls not only single-cell dy- 
namics but also the strength and time scales of spike correlations in cell 
groups receiving common inputs. We show how this relationship can be un- 
derstood. As a membrane current (here, g^) is adjusted, firing rate - current 
curves progressively transition (here, from Type I to Type II). At the same 
time, there is a transition in periodic orbit types that neural trajectories 
visit (here, ranging from orbits "near" a fixed point to relatively "isolated" 
orbits (Rush and Rinzel, 1995)). In turn, this produces a steady progression 
of spike-triggered averages — and hence the filters that neurons apply to 
shared input signals (here, from primarily integrating to primarily differen- 
tiating modes (Mato and Samengo, 2008), cf. ( Agiiera y Areas et al., 200~3| )). 



Basic formulas can then be used to translate these filtering properties into 



predictions for correlated spiking in neural pairs and populations (Ostojic 



et al., 2009) as well as the downstream impact of this cooperative activity. 
We anticipate that this approach will bear fruit in studies of the collective 
activity of a wide variety of neuron types. 



Relationship with prior work 

A number of prior studies have considered the problem of how spike generat- 
ing dynamics affect the transfer of incoming current correlations into outgoing 



spike correlations (Binder and Powers, 2001[ de la Rocha et al., 2007| Shea 



Brown et al., 2008] Rosenbaum and Josic, 20lTj Tchumatchenko et al., 2010] 



Marella and Ermentrout, 2008 


Vilela and Lindner, 2009 


Barreiro et al., 


2010 


Ostojic et al., 2009 |Tchumatchenko et al., 2010 


Hong and De Schut- 



ter, 2008). In particular, flde la Rocha et al., 2007; Shea-Brown et al., 2008 



Rosenbaum and Josic, 2011 ) show that leaky integrate-and-fire (LIF) neurons 
can transfer up to 100% of current correlations into spike count correlations. 
The level transferred increases with the firing rate at which single neurons 
are operating, and the time scale T. These findings are simpler to state 
compared with the present results for conductance-based neuron models, for 
which 100% correlation transfer is never obtained, and trends with T differ 
depending on g&. 



Other works (Hong and De Schutter, 2008 Shea-Brown et al., 2008 
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Vilela and Lindner, 2009) investigate correlation transfer in more complex 



spiking models. In particular, Shea-Brown et al. (2008) and Vilela and Lind- 



ner (2009) explore the full parameter space of input currents for the quadratic 

- arguably that with the next level of complexity 
These authors find similar trends in correlation 



integrate and fire model 
beyond the LIF model, 
transfer as the neurons' operating points change, but a limitation to 66% 



rather than 100% in correlation transfer. Meanwhile, Hong and De Schut- 
ter (2008)] show complex dependencies on neural operating point for the 
Hodgkin-Huxley model. Taken together, these studies suggested that corre- 
lation transfer depends on spike generating dynamics in a rich and diverse 
ways. 

This opened the door to a broader study, but exploring correlation trans- 
fer for the full space of possible spike generating dynamics in neural models is 
a daunting task. The axis that spans from Type I to Type II excitability pro- 
vides a natural focus. This has been explored using sinusoidal "normal form," 
phase-reduced models (Marella and Ermentrout, 2008 Galan et al., 2007 



Barreiro et al., 2010 Abouzeid and Ermentrout, 2011). These studies used 



simulations in the superthreshold regime, together with analysis in the limits 
of very short or very long time scales T, to show the same trend in correla- 
tion transfer over short vs. long T that we find here for conductance-based 
models. A greater frequency of instantaneous (small T) spikes for Type II 



vs. Type I neurons was predicted using these simplified models (Marella and 



Ermentrout, 2008; Galan et al., 2007); later, Barreiro et al. (2010) predicted 



the switch in relative correlation transfer efficiency from Type II to Type I 
models as T increases. 

The present contribution is to test the resulting predictions using bio- 
physical, conductance-based models valid in a wider range of firing regimes, 
to explain the origins of variable correlation transfer via filtering properties 
of Type I vs. Type II cells, and to demonstrate the impact on downstream 
neurons. 



Scope, limitations, and open questions 

The circuit model that we have studied, as illustrated in Fig. [TJ is limited 
to a single, idealized feature of feedforward connectivity: overlapping in- 
puts to multiple recipient cells. More realistic architecture could include 
delays in incoming inhibitory vs. excitatory inputs (Gabernet et al., 2005). 



Interactions of shared-input circuitry with recurrent connectivity also pose 
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important questions (Ly and Ermentrout, 2009). This is especially so given 
the distinct properties of Type I vs. Type II cells in synchronization due to 



reciprocal coupling (Rinzel and Ermentrout, 1998; Ermentrout and Terman. 



2010). 



Other aspects of our biophysical and circuit dynamics are also idealized. 
For one, individual input currents fluctuated on arbitrarily fast time scales 
(i.e., as white noise processes). Relaxing this would be an interesting ex- 
tension. While prior studies (de la Rocha et al., 2007) suggest that trends 
will persist for inputs with fast (but finite) time scales, new effects could 
arise for slower-time scale inputs representative of slower synapses or even 
network-level oscillations. Another addition would be for inputs to arrive via 
excitatory and inhibitory conductances, rather than currents; while previous 
studies with integrate and fire cells (de la Rocha et al., 2007) have found 
that this yields qualitatively similar results, there could be interesting inter- 
actions with underling filtering properties in biophysical models. The same 
holds true for inputs that arrive at dendrites in multi-compartment models. 

Likewise, the circuitry of Fig. [To|a) that we used to investigate the im- 
pact of correlated spiking on downstream neurons was highly idealized. An 
especially appealing extension would be to note that inhibitory and excita- 
tory neurons often have different excitability types. Thus, downstream cells 
could receive input from both excitatory Type I and inhibitory Type II pop- 
ulations. Our results suggest that sensitivity to excitatory vs. inhibitory 
afferents would vary with membrane time constants downstream, possibly 
amplifying the modulatory effects identified here. 
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